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ABSTRACT 

In the framework of the concordance cosmological model the first-order scalar and 
vector perturbations of the homogeneous background are derived in the weak gravi¬ 
tational field limit without any supplementary approximations. The sources of these 
perturbations (inhomogeneities) are presented in the discrete form of a system of sep¬ 
arate point-like gravitating masses. The found expressions for the metric corrections 
are valid at all (sub-horizon and super-horizon) scales and converge at all points ex¬ 
cept at locations of the sources. The average values of these metric corrections are zero 
(thus, first-order backreaction effects are absent). Both the Minkowski background limit 
and the Newtonian cosmological approximation are reached under certain well-defined 
conditions. An important feature of the velocity-independent part of the scalar pertur¬ 
bation is revealed: up to an additive constant this part represents a sum of Yukawa 
potentials produced by inhomogeneities with the same finite time-dependent Yukawa 
interaction range. The suggested connection between this range and the homogeneity 
scale is briefly discussed along with other possible physical implications. 

Subject headings: cosmological parameters — cosmology: theory — dark energy — dark 
matter — gravitation — large-scale structure of universe 


1. INTRODUCTION 


The concordance cosm o logical model fits we ll the contemporary observations (jHinshaw et al 


201.11 : lAde et al.l 120141 . l201,T lAubourg et al.ll201,^ ). This model assumes that the Universe is filled 
with dominant portions of cold dark matter (CDM) and dark energy, represented by the cosmo¬ 
logical constant and assuring the acceleration of the global expansion, as well as a comparatively 
small portion of standard baryonic matter and a negligible portion of radiation. According to the 
cosmological principle, on large enough scales the Universe is treated as being homogeneous and 
isotropic, so the corresponding background Friedmann-Lemaitre-Robertson-Walker (FLRW) met¬ 
ric is appropriate for its description. On the contrary, on sufficiently small scales the Universe is 
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highly inhomogeneous (separate galaxies, galaxy groups and clusters stare us in the face). The 
thorough theoretical study of the structure formation starting from primordial fluctuations at the 
earliest evolution stages with the subsequent comparison of the predictions with the cosmic mi¬ 
crowave background and other observational data may be recognized as one of the major subjects 
of modern cosmology. 


The structure growth is usually investig ated by m e ans o f two main disti n ct approaches, namely, 
the re lativistic perturbation theory (see, e.g.. lBardeenl (|l98Cll i: lDurreii (j2008l i : iGorbunov &: Rubakov 
(|201l|)) a nd A^-body sim u lations generally b a sed on the Newtonian cosmo logical approximation 
(see, e.g., ISoringell (j2005l i: lOolag et al.l (|2008l i: IChisari fc Zaldarriagal ((20111)). Roughly speaking, 
the area of the first approach may be characterized by the keywords “early Universe; linearity; large 
scales”, while the area of the second one may be defined by the keywords of the opposite meaning: 
“late Universe; nonlinearity; small scales”. Both approaches make great progress in describing the 
inhomogeneous world within the limits of their applicability. Nevertheless, the linear cosmological 
perturbations theory certainly fails in describing nonlinear dynamics at small distances, while 
Newtonia n simulations do not take into account rela t ivistic effects becoming n on-negligible at large 


distances (Green fc Wald 


2012 


Adamek et al 


2013 


20141 : iMilillo et al.l l2015l i . In this connection 


in the latter case certai n effort is required for extracting relativistic features from the large-sca le 


Newtonian description (jGhisari fc Zaldarriaga 


2011 


Fidler et al.ll2015l : iHahn Sz ParaniaDell2016l ). 


Until now there was no developed unified scheme, which would be valid for arbitrary (sub¬ 
horizon and super-horizon) scales and treat the non-uniform matter density in the non-perturbative 
way, thereby incorporating its linear and nonlinear deviations from the average. This acute problem 
is addressed and successfully resolved in the present paper by construction of such a self-consistent 
and indispensable scheme, which promises to be very useful in the precision cosmology era. 

A couple of previous similar attempts deserve s mentioning. First, the g eneralization of the well- 


known nonrelativistic post-Minkowski formalism ((Landau fc Lifshitz 


2000) to the cosmological case 


in the form of the relativistic post-Friedmann formalism, which would be valid on all scales and 
include t he full nonliii e aritv of Newtonian gravity at small distances, has been made in the recent 
paper bv IMilillo et al.l ((2010), but the authors resorted to expansion of the metric in powers of the 
parameter 1/c (the inverse speed of light). Second, the formalism for relativistic A"-b ody simulations 


i n the weak field regime, suitable for cosmological applications, has been developed bv lAdamek et al 


((2010), but the authors gave different orders of smallness to the metric corrections and to their 


spati al derivatives (a similar “dictionary” can be found in ((Green fc Waldl (20121 : (Adamek et al 

2o3)). 


The current paper also relies on the weak gravitational field regime: deviations of the metric 
coefficients from their background (average) values are considered as first-order quantities, while 
the second order is completely disregarded. However, there are no additional assumptions: in 
the spirit of relativity, spatial and temporal derivatives are treated on an equal footing, and n o 
dictionary giving them different orders of smallness is used (in contrast to (Green &: Waldl (|2012l i ; 
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Adamek et al.1 (j2013. l2014l'lk expansion into series with respect to the ratio 1/c is not used as 


well (in contrast to iMilillo et alJ (j2015l ii: there is no artificial mixing of first- and second-order 
contributions; the sub- or super-horizon regions are not singled out, and the derived formulas for 
the metric corrections are suitable at all scales. 


_ The desired formalism is elaborated below within discrete cosmo logy (jEingorn fc Zhukll2f)12l . 

2014 ; Eingorn et al.ll2013 ; Gibbons fc Ellis 2014 ; Ellis h GibbonsIboiS ), based on the well-grounded 
idea of presenting nonrelativistic ma tter in the form of separa te point-like particles with the corre¬ 
sponding energy-momentum tensor ([Landau &: Lifshita 120001 1. At sub-horizon scales, implementa¬ 
tion of this idea leads to nonrelativisti c gravitational potentials and Newtonian equations of motion 
against the homogeneous background (|Peebleslll980l L which are commonly used in modern A^-body 
simulations. Now the explicit expressions for potentials, applicable at super-horizon scales as well, 
will be made available. 


The paper is structured in the following way. Section 2 is entirely devoted to solving the 
linearized Einstein equations for the first-order scalar and vector cosmological perturbations of the 
homogeneous background. In Section 3, the arresting attention properties of the derived solutions, 
including their asymptotic behavior, are analyzed and their role in addressing different related 
physical challenges is indicated. The main results are summarized laconically in the Gonclusion. 


2. DISCRETE PICTURE OF COSMOLOGICAL PERTURBATIONS 

2.1. Equations 


The unperturbed FLRW metric, describing the Universe being homogeneous and isotropic on 
the average, reads: 

ds^ = {drf' — , a, /3 = 1, 2,3 , (2.1) 

where a{r]) is the scale factor; r] is the conformal time; x“, a = 1,2,3, stand for the comoving 
coordinates, and it is supposed for simplicity that the spatial curvature is zero (the generalization 
to the case of non-flat spatial geometry is briefly analyzed below). The corresponding Friedmann 
equations in the framework of the pure AGDM model read: 




= K£ + K 


( 2 . 2 ) 


and 


m' + 


= A. 


(2.3) 


where Ti. = a'/a = {da/dr])/a] the prime denotes the derivative with respect to rj; k = SttGn/c^ 
(c is the speed of light and Gjsf is the Newtonian gravitational constant); e represents the energy 
density of the nonrelativistic pressureless matter; the overline indicates the average value, and A 
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is the cosmological constant. Domination of cold matter and A is at the center of attention, so 
contributions of radiation, relativistic cosmic neutrinos or any warm component are negligible. 


Fo dowing the analysis of t he fir st-order cosmological perturbations by iBardeenl (jlQSOl ) ; iDurrer 


(120081 ) : lGorbunov fc RubakovI (j201ll ). let us fix the Poisson gauge and consider the respective metric 


ds^ = m 


(1 + 2 <^) drf' + 2Badx°'dr] — (1 — 2 <^) Sapdx^dx^ 


(2.4) 


where the function <^(ry,r) and the spatial vector B(r/, r) = (Bi, B 2 , B^) describe the scalar and 
vector perturbations, respectively. It is assumed that there is no anisotropic stress acting as a source 
of the difference between perturbations of the metric coefficients goo and gap, a = (3. Therefore, 
both these metric corrections are equated to the same expression 2o^<^ from the very beginning. 
Tensor perturbations are not taken into account because their source is also of the order beyond 
the adopted accuracy. The first-order tensor perturbations are associated with gravitational waves, 
freely propagating against the FLRW background. Their propagation is governe d by the well- 
know n equation uncoupled from the equations for 4> and B (see, e.g., Eq. (11) in (|Noh &: Hwang 
200a)). Here attention is concentrated solely on the perturbations with n on-negligible sources , 
and c osmological gravitational waves are not investigated. Similarly to, e.g.. lAdamek et al.l (l2013l . 
2 OI 4 I I. it is demanded that 

(2.5) 


VB = = 0. 


Then the Einstein equations 




G^ = KTt + A6f, i,k = 0,l,2,3, 


( 2 . 6 ) 


where and denote the mixed components of the Einstein tensor and the matter energy- 
momentum tensor, respectively, take the following form after linearization: 


G^ = kTI^ + A => A<l>+ n^) =-na^dT^, 


GO = kT° 




(2.7) 

( 2 . 8 ) 


G^ = + A6^ 


4>" + 3 n ^' + { 2 H ' + 7^2^ $ = 0 , 


(2.9) 


dBa ^ dB, 


dxd dx°‘ 


+ 2'H 


I dxf^ ' 


dBa ^ dB 


= 0 . 


Here A = 5"^ 


( 2 . 10 ) 


92 

■ stands for the Laplace operator in the comoving coordinates; +STk, 


dx°‘dx^ 

where Tg = e is the only nonzero average mixed component. In the spirit of the particle-particle 
method of A-body simulations, the matter constituent of the Universe may be presented in the 




form of separate point-like massive particles. Then the deviations 5Tk from the average values Tj 
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can be easily determined with the help of the well-known general expression for the matter energy- 


momentum tensor con t ravar iant components ([Landau &: Lifshital200Cll : IChisari &: Zaldarriagall2011 
Eingorn &: Zhuk 2012 . 2014 1: 


rj,ik ^ dx\, dxj dr] ^ 


dr] dr] dsr. 


( 2 . 11 ) 


where g = det{gik)- This expression corresponds to a system of gravitating masses with the 
comoving radius-vectors r„(7/) and 4-velocities = dxl^jdsn, i = 0,1, 2,3. Their rest mass density 
p{g, r) in the comoving coordinates reads: 


P = '^mn5{r -Tn) = '^Pn, Pn ^ mnd(r - Fn) ■ 


( 2 . 12 ) 


Introducing the comoving peculiar velocities 1)“ = dx°/dg, a = 1,2,3, and treating them as 
importing the first order of smallness in the right-hand side (rhs) of the linearized Einstein equations 
(1231) and (|2.8p . one finds out that within the adopted accuracy 


5Tq =Tq — Tq — -^5p -\ 5p — p — p, 




2>pc^ 


(2.13) 


and 


c2 ^ 


<5r“ = —3 ^ m„<5(r - r„)u“ — 


— 2 9 — 2 

^ pC^ ^ , P(^ r-. 

~ j Pn'^n ~^Bq 


(2.14) 


respectively. In addition, 6Ti = 0. In particular, in full agreement with, e.g., Adamek et al. ( 2013 1. 
the first-order anisotropic stress is considered as vanishing for nonrelativistic matter. Consequently, 
the rhs of both Eqs. ()2.9I) and (12.101) are zero. In the formula (12.131) the arising term ~ p^ is replaced 
by the term ~ since the product 6p^ imports the second order of smallne ss and therefore should 
be dr opped (the inequality \5p\ S> certainly holds true at all scales (jChisari fc Zaldarriaga 

2 OIII II. The same reasoning applies to the term ~ pBa arising in the formula (|2.14p . The average 
rest mass density p is related to the average energy density e by means of an evident equality: 
e = 'pr? jc?. 


It is crucial that t hroug h out the presen t paper, sim ilarly to IChisari &: Zaldarriagal (|201lll : 
Eingorn fc Zhuk ( 2012 . 2014 1: Adamek et al. ( 2013 . 2014 1. the rest mass density p is treated in 
the non-perturbative way: fulfilment of the inequality \5p\ <C p is not demanded. Eor instance, 
the intragalactic medium and dark matter halos are characterized by the values of p being much 
higher than p. Thus, nonlinearity with respect to the deviation of p from its average value p at 
small scales is fully taken into consideration here. The sole requirement lies in the following: this 
deviation 5p as a source of the metric correction must secure smallness of this correction, i.e. 
the inequality |4>| <C 1. As regards the Einstein tensor components, their nonlinear deviation from 
the corresponding background values is also unforbidden. For instance, the term ~ A<h in the 
expression for Gg predominates at small scales not only over the other terms in this expression, but 
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also over the background value /a? of Gq. Nevertheless, as distinct from I Adamek et alJ (2013, 
201J), those terms which are nonlinear with respect to the metric corrections are neglected in the 
expressions for the Einstein tensor components. This concerns, e.g., the product <1>A<1> which is 
neglected in comparison with Ad>, in complete agreement with the initial well-grounded assumption 
|d>| < 1 . 


Substitution of the expressions ()2.13p and (I2.14p into Eqs. (|2.7p and (12.81) . respectively, gives 

..2 


A$ - m{^' + n^) - = '^5p , 

la la 


-AB + V($ + m) - —B = -^ E = -^ E 


Pn^n \ 


where ^nip) = dxn/dr] = With the help of the continuity equation 

Pn + '^iPn^n) = 0 , 


(2.15) 

(2.16) 

(2.17) 


which is satisfied identically for any n-th particle, it is not difficult to split the vector ^ PnVn 

n 

into its grad- and curl- parts: 


Pn^n = '^^ + Pn^n - VH , 


where 


-V 

A'TT * ^ 


r - r„ 


m. 


dvr ^ Ir-FnP 


Really, this function satisfies the Poisson equation 

AS = V y] PnVn = - y] Pn > 

n n 

and this fact can be easily checked alternatively with the help of the Fourier transform: 

-A:=|kl, 


where 


1 ( 77 , k)= / S(r 7 , r) exp(—ikr)fir , 


(2.18) 

(2.19) 

( 2 . 20 ) 

( 2 . 21 ) 

( 2 . 22 ) 


P„(t 7 , k) = / p„(t 7 , r) exp(-ikr)dr = mn (5(r - r„) exp(-ikr)dr = rUn exp(-Akr„). (2.23) 


In other words, one can demonstrate that the function ^(ry, k) derived from (12.2211 after the substi¬ 
tution of (I2.19P satishes Eq. p2.2ip and, consequently, reads: 


p y~]mn(kv„)exp(-Akr„). 


( 2 . 24 ) 
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Taking into account (|2.18p . from (I2.16P one gets 





and 



— 2 2 

npc g nc 

2a 2a 


Substitution of (I2.25P into p2.15p gives 



A$ - 


3Kpc^ 

2a 


$ =- 6p 

2a ^ 


3acc^77„ 

2a 


(2.25) 


(2.26) 


(2.27) 


Thus, Eqs. (I2.26P and p2.27p are derived for the vector and scalar perturbations, respectively. Below 
their solutions are found and the fulfilment of Eqs. (12.5p . (12.9p . (I2.10p and (I2.25j) is verified. 


2.2. Solutions 

In the Fourier space Eq. (|2.26l) takes the following form: 

--B-B =- y Pn^n — 

4 2a 2a 

Substituting ()2.23p and (12.241) into p2.28p . one immediately obtains 

B = 


2ac^ /.9 2 k/0c^^ ^ 


k^ + 


(2.28) 


rrin exp(-zkr„) ( v„ - 1 . 


(2.29) 


The condition (j2.5p is evidently satisfied since kB = 0. It is also not difficult to verify that 
Eq. (I2.10p is fulfilled within the adopted accuracy. In fact, it is enough to show that 


B' + 277B = 0 . 

For this purpose let us write down the spacetime interval for the n-th particle 


dSf] — a 


nl/2 


1 + 2$ + 2Bau“ - (1 - 2$)(5„/3u“u(^ dp 


(2.30) 


(2.31) 


and the corresponding Lagrange equations of motion 


[a(B|r=r„ - V„)]' = oV$|r=r„, 


(2.32) 


where the contributions of the considered particle itself to B and $ are excluded as usual, so 
divergences are absent (in other words, the particle moves in the gravitational field produced by 
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the other particles). Since Eq. (I2.32p and its consequences will be used exclusively in linearized 
Einstein equations (e.g., Eq. (12.101) 1. all terms being nonlinear with respect to B and <h have 
been dropped in order to avoid exceeding the adopted accuracy. Multiplying ()2.32p by pn and 
summing up, one gets 

P(oB)'- = apV^. (2.33) 

n 

Further, in the terms containing B and $ the rest mass density p should be replaced by its average 
value p as discussed before. Consequently, 

'^Pn{avn)' = -apV$ + p(aB)'. (2.34) 

n 

Of course, in order to study dynamics of the A^-body system one should include the removed 
nonlinear terms —a5p\/^ and 5/?(aB)h However, for analyzing the Einstein equations for the first- 
order cosmological perturbations it is apparently enough to keep only linear terms in the equations 
of motion. In the Fourier space Eq. (j2.34p reads: 

^/5n(av„)'= ^w-nexp(-ikr„)(avn)' = -op • zkl>p(aB)'. (2.35) 


Now, expressing B' from (I2.29P with the help of (j2.35l) and substituting the result into (I2.30p . one 
arrives at the identity. 

Finally, the vector perturbation B can be determined by multiplying (12.291) by exp(ikr)/(27r)^ 
and integrating over k. Cumbersome, but straightforward calculation gives 


B = —Y" 

Svro ^ 


rrin^n {3 + 2V3qn + ‘iql)exp{-2qn/V3)-3 


r - r. 


qI 




m.n[vn(r - r„)] 


r - 


(r - rn) 


9 - (9 6\/3gn + ^n) exp(-2q'„/V3) 


(2.36) 


where the following convenient spatial vector is introduced: 


qn(p,r) 



I’n); Qn — |qn 


(2.37) 


Let us now switch over to Eq. (I2.27p . In the Fourier space this equation takes the following 


form: 


, 2 i 3 kpc"‘ 
- k ^ - 


2a 


i =—y 

2a ^ 


Pn 




2a 


2a 


(2.38) 


where the well-known presentation (27r)^(5(k) = f exp(—ikr)dr is taken into account. Substituting 
(j2.23p and (|2.24l) into (|2.38p . one immediately obtains 


KC 




3kp& 

2a 


2 \ -1 


nin exp(—ikr„) ( 1 -|- SiTi 






- p{2Trf6{k) 


(2.39) 
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With the help of (|2.35p one can show that both Eqs. (12.91) and (|2.25p are fulhlled within the 
adopted accuracy. The scalar perturbation $ itself can be determined by multiplying (j2.39l) by 
exp(ikr)/(27r)^ and integrating over k. This Fourier inversion gives 


$ = 


+ 


1 

3 



n 


mn 

r - Tn 


exp(-g„) 


3^.^ mn[vn{r - r^)] _ 1 - (1 + gn)exp(-g„) 
Svra 4^ |r-r„| 


(2.40) 


Thus, the explicit analytical expressions (12.361) and (|2.40l) for the first-order vector and scalar 
cosmological perturbations, respectively, are determined for the first time. Let us accentuate the 
irrefutable fact that the dictionary-based approach compels us to say goodbye to all hopes of finding 
analytical solutions. Really, let us momentarily return to Eq. (|2.27l) . In (jAdamek et al.l l2013l l 
the order 0(e) is assigned to $ while every spatial derivative is treated as importing the order 
0(e“^/^), where e is some small parameter. Consequently, A<I> ~ 0(1) while [3Kpc^/(2a)]<h 0(e). 

Concerning the rhs of Eq. (I2.27p . following the same assignment, 6p ~ 0(1), but S ~ 0(e) (see 
Eq. (l2.25p L Thus, in the 0(l)-approximation the second terms in both sides of Eq. (12.271) are 
missing, and this equation reduces to the Poisson equation of Newtonian gravity, being inapplicable 
for large enough distances. lAdamek et al.l ([2013) overcome this difficulty by including the terms of 
the orders 0(1) and 0(e) simultaneously in the same equations (i.e. by mixing their own 0(1)- and 
O(e)-approximations). Then the term [SKpc^ / (2a)]4> is reinstated, but so do the nonlinear terms like 
(1A4> ~ 0(e). This substantially complicates the problem leaving room for numerical computation 
only. As shown in the current paper, it is mathematically logical to keep both considered terms 
A<I> and [3Kpc^/(2a)]4> (linear with respect to 4>), dropping all nonlinear ones such as <I>A4>. The 
resulting Helmholtz equation (|2.27p is not restricted to sufficiently small distances and actually 
covers the whole space. Therefore, the desired formalism for cosmological iV-body simulations is 
developed here shedding hardly any blood: linearized Einstein equations are not complicated by 
extra nonlinear terms and admit exact analytical solutions describing contributions of every single 
massive particle to the total inhomogeneous gravitational field. 


In what follows the noteworthy features and advantages of the expressions (|2.36p and (|2.40p 
are briefly discussed. 


3. MENU OF PROPERTIES, BENEFITS, AND BONUSES 
3.1. Minkowski Background Limit 


After obtaining the solutions (12.361) and (|2.40p of the system of Eqs. (|2.7p - (l2.10p for the hrst- 
order cosmological perturbations, it is absolutely necessary to study their asymptotic behavior. 
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First of all, let us consider the Minkowski background limit: 

a —)■ const ^ )-0; >-0 = 

Then instead of ()2.36p and (I2.40p one has, respectively, 

m„.[v„(r - r„)] 


9n 0 . 


(3.1) 


B 


KC 

dvra 

n 

2 c 2 


rrinVr, 


r - 


+ 


r - 


(r - Tr, 


E 


rrin 


R — Rri 


4v„ + 


4[v„(R-R„)] R-R„ 


|R — Rri 


|R — Rri 


where the physical radius-vectors R = ar, R„ = ar„ are introduced, and 

KC^ ^ nin Gn 

$ -i,-\ 

Svra ^ 


r - Tr, 


Gn 

'TF 


nir. 


R — Rri 


(3.2) 


(3.3) 


where the constant 1/3 has been dropped since it originates in (I2.27P exclusively from the terms 
containing p. Let us compare the asymptotic expressions (13.21) and (13.31) with the corresponding 
corrections to Min k owski spacetime metric coefficients from the paragraph 106 of the textbook by 
Landau fc Lifshita (120001 ). This paragraph is devoted particularly to the weak gravitational field 
generated by a system of nonrelativistic point-like particles perturbing flat spacetime geometry. 
Multiplying (13.3p by 2, one arrives at the resul t which exactly coincides with the first term in the 
expression (106.13) in (jLandau fc Lifshital2000l l for the metric correction h^Q, as it certainly should 
be. In addition, the term containing v„ disappears from (I2.40p in the considered limit (in view of 
the factor T-L), and at the same time there is no term that is linear with respect to in (106.13) 
([Landau &: Lifshital2000l L This fact also serves as confirmation of coincidence. 


As regards (13.2p . the only difference between this expression and (106.15) in ([Landau fc Lifshitz 


2000|) is that the integers 4; 4 are replaced by 7; 1 respectively (one should keep in mind that the 
comoving peculiar velocities defined with respect to the conformal time rj as dvn/dr] are related 
to those defined with respect to the synchronous time t as Vn ^ dr^/dt by means of an evident 
equality: cdt = adrj ^ v„ = avn/c). Apparently, this difference in integers rep resents none 
other than a result of different gauge conditions here and in ([Landau &: Lifshitzjl2000lL Indeed, th e 
condition (12.5p is not demanded and, of co urse, does not hold true in ([Landau Sz Lifshita l2000lL 
Correspondence between (13.2p and (106.15) ([Landau fc Lifshital200Clf lies in the fact that the sum 
of these integers is the same: 4-|-4 = 7-1-1. One can show that it equals 8 for the other appropriate 
gauge choices as well. T herefore, (|3.2[) exactly coincides with the purely vector part of (106.15) 
([Landau &: Lifshitjl200ril ). as one can easily see by finding curl of both these expressions. 


3.2. Newtonian Approximation and Homogeneity Scale 

Now let us switch over to the Newtonian cosmological approximation: <C 1, i.e. |r — r„| <C 

y^2a/{3Kpc^), and peculiar motion as a source of the gravitational field is completely ignored 
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(jChisari &: Zaldarriagall201lh . so the summands directly proportional to the velocities are omit¬ 
ted. Then only the scalar perturbation <1> survives in the same form (13.3p . where the constant 
1/3 has been dropped for the other reason: only the gravitational potential gradient enters into 
equations of motion describing dynamics of the considered system of gravitating masses. These 
equations for any j-th particle follow directly from (I2.3ip and take the form 


R-i —Rj = —gn 'y 

a ^ 


rrir, 


(R-i — Rr, 


|Rj ~ Rri 


(3.4) 


in the physical coordinates 0 = 1,2,3, being in ac c ordance with t h e correspondin g 

equations in t h e papers by Soringel (2005); Dolag et al. (2008); Labini ( 2013 ): Warren ( 2013 1: 
Eingorn (2014); Ellis h GibbonsI (2015) devoted to cosmological simulations. Here dots denote the 
derivatives with respect to t. 

Let us consider two important questions. First, what are the applicability bounds for the 
above-mentioned inequality, which may be rewritten in the form |R — R„ I <C 20^/{‘inpc^)l In 
order to answer, one should simply calculate the rhs of this inequality: 




I 2a? 
2>KpC^ 


2c^ 


9H^nM 


ao 


KpC 

3^’ 


(3.5) 


where oq and Hq are the current values of the scale facto r a and the H ubble parameter H = a/a = 
{da/dt)/a = cH/a, respectively. According to Ade et al. ( 2014 . 2015 ). Hq k, 68kms“^Mpc“^ and 
Q-m ~ 0.31. Therefore, the current value of A is Ao ~ 3700 Mpc « 12 Gly. It is very in teresting that 


this Yukawa interaction range and the sizes of the largest known cosmic structures (j Clowes et al 


2013; 


Horvath et al.ll2014l : iBalazs et al.ll2015l ) are of the same order, thereby hinting at the oppor¬ 


tunity to resolve the formidable challenge lying in the fact that their siz es essentially exce ed the 
previously reported epoch-independent scale of homogeneity ~ 370 Mpc (|Yadav et al.l 120101) . The 
authors arrived at this underestimate by comparing the deviation of the fractal dimension, char¬ 
acterizing the distribution of matter, from 3 (dimensionality of space) to its statistical dispersion. 
Along with fractal analysis, their approach relies on the weak clustering limit and cosmological 
simulations driving 512^ particles in a cube with the edge ~ 1.5 Gpc. Incidentally, this edge is 
less than half Aq, and from the very beginning of such a volume-restricted simulation it is diffi¬ 
cult to expect any definite and reliable indications of structuring in bigger volumes. Now, if one 
associates the scale of homogeneity with A instead, then the cosmological principle, asserting that 
the Universe is homogeneous and isotropic when viewed at a sufficiently large scale, is saved and 
reinstated when this typical averaging scale is greater than A. The proposed association does not 
mean that the homogeneity scale is equated exactly to A but rather describes A as an approximate 
upper bound to the cosmic structure size, and the homogeneity scale as a distance exceeding A 
in a few tim es while remaini ng of the same order. It is remarkable that this reasoning is actually 
confirmed bv iLi fc: LinI (j2015l L The authors dehned the scale of homogeneity as a distance at which 
the correlation dimension is within 1 % of 3 (and, consequently, equals 2.97) and fixed an upper 
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bound to such a distance ~ SAq. The dependence A ~ is noteworthy as well: the earlier 
the evolution stage, the smaller the scale of homogeneity. Natnrally, this is closely related to the 
hierarchical clustering process. 

The second important question is: what are the applicability bounds for peculiar motion 
ignoring? In order to answer, one can consider the ratio of the third term in (12.401) to the second 
one. For a single gravitating mass mi momentarily located at the origin of coordinates (ri = 0) with 
the velocity vi collinear to r (for ensuring the maximum value of the scalar product vir = v^r, where 
hi = |vi|) this ratio amounts (up to a sign) to 2>'Hvir/2 = 3HaviR/ (2c^), where vi = [vil = cvi/a, 
R = |R| = or and gi <C 1 as before. Actually the product avi is none other than the absolute 
value of the particle’s physical peculiar velocity. For example, with the help of the today’s typical 
values (250 4- 500)kms“^ and the inequality R <C 3700 Mpc one finds that the considered ratio is 
much less than (1 2) x 10~^. Exactly the same estimate can be made for the ratio of derivatives 

of the considered terms with respect to r. This means that at the scales under consideration, the 
gravitational force originating from the second term in the gravitational potential (j2.40jl . which 
does not contain particle velocities, is much stronger than that coming from the third one, which 
contains them. 

Thus, the Newtonian cosmological approximation may be used when |R—R„| <C A. Otherwise, 
at the scales comparable or greater than A, one should use the complete expressions for the metric 
corrections obtained in the previous section. In particular, the derived Yukawa-type potentials 
should be used instead of the Newtonian ones in order to study formation and evolution of the 
largest structures in the Universe. It is necessary to understand that the elaborated formalism 
results in Newtonian behavior of the considered physical system at sufficiently small distances 
without any relativistic corrections. Hence, the accuracy of the developed theory is limited in this 
region by the standard Newtonian approach. However, the predicted Yukawa behavior at greater 
distances may be considered a relativistic effect since it follows directly from Einstein equations of 
General Relativity. 

It should be emphasized that despite the presence of those terms in (I2.36|) and (I2.40p . which 
do not contain exponential functions, the influence of any particle on the motion of its neighbours 
does drop exponentially when the distance increases. Really, with the help of (I2.30p the equations 
of motion (I2.32p may be rewritten in the form 

(av„)' = -a(V4>|r=r„+77B|r=r„) , (3.6) 

and this pecnliar acceleration of a given particle, caused by all other gravitating masses, contains 
solely terms with exponential functions. Indeed, the direct substitution of ()2.36l) and ()2.40p into the 
rhs of (13.6p demonstrates that all terms without exponential functions exactly cancel each other. 
This fact confirms the revealed Yukawa nature of universal gravitation. In addition, it indicates 
that for sufficiently small values of a and nonzero separation distances between particles they al¬ 
most do not interact gravitationally (all terms containing exponential functions may be dropped 
under such conditions), so the system behaves as a perfect gas undergoing the global expansion. It 
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is also interesting that the physical screening length \/3A/2 from ()2.36p is less than the counterpart 
A from (|2.40p meaning that vector modes diminish with distance faster than scalar modes. The 
equations of motion (13.61) are ready to be used in a new generation of cosmological simulation codes 
(see, however, the discussion of coordinate transformations below, indicating possibilities of rein¬ 
terpretation of Newtonian simulations from a relativistic perspective). It would be quite reasonable 
to confront the outputs of relativistic simulations with those of various Newtonian predecessors, 
thereby discriminating between them and the proposed Yukawa modification, especially regarding 
predictions of peculiarities of hugest gravitationally bound objects in the Universe. 

One more important detail consists in the fact that A does not coincid e with the Hubble 
radius c/H, in contrast to the Yukawa interaction range proposed bv ISignord (|2005l l in order to 
limit gravitational effects of a particle outside its causal sphere. Really, in terms of the Hubble 
parameter H and the deceleration parameter q = —a/ 


1 3^2^^ ^ 3H 

A2 “ + “ "lY’ 


A = 


\/3(l -I- q) 


H 


(3.7) 


following directly from (13.Sp and the Friedmann equations (12.2p and (|2.3I) . which may be rewritten 
in the form 


3772 


K/9C 


+ A 


and 


o ^ 3772 ^ 277 ^ 

^(1 - 2q) = - 5 - = A, 


(3.8) 


(3.9) 


respectively. One obtains from (|3.7I) that A = c/77 in the unique moment of time when q = —2/3. 
According to Eqs. ()3.8p and (I3.9D. 2A/7 = Kp c ^/a^ ai this moment, or 2Q\/7 = 71 ^( 00 / 0 )^) 
where Ua = Ac^/ (377q) « 0.69 ( Ade et al. 2014 . 2015 1. Hence, A = c/H in the near future when 
a/oo ~ 1.16. Before this moment A < c/H, while afterwards the opposite inequality takes place. 


Likewise A does not coincide with a shielding length introduced bv lHahn &: Paraniaoel (|2016l l. 
The authors resorted to the dominant growing mode in the framework of the linear relativistic 
perturbation theory (see their Eq. (15), which is actually a predetermined approximate solution 
but, nevertheless, serves as an assumed starting point) and presented <I> in the standard form of 
a product of a function of time and a function of spatial coordinates. This allowed expressing 
377 (4>' -|- 774>) as where I is a certain time-dependent parameter, and then, after substitu¬ 

tion into the linearized Einstein equation Gq = kT// + A, declaring / to be a shielding length. It 
should be mentioned that the s ame s hielding mechanism may be also discerned in the preceding 
paper by lEingorn fc BrilenkovI (j2ni5l ) , where continuous matter sources are in the attention fo¬ 
cus instead of discrete ones investigated here. In this connection, it makes sense to confront in 
brief the approaches bv iHahn &: Paraniaoel (|2016l l and lEingorn &: BrilenkovI (|2015l l. Eirst, at the 
same level of linear e nergy -momentum fluctuations, the velocity-dependent term introduced by 
Eingorn &; BrilenkovI (j2015l l in the equation for (see, e.g., their Eq. (16)) can be also easily re¬ 
duced to for the considered growing mode. Of course, this is a foreseeable coincidence because 
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the mentioned velocity-dependent term coincides exactly with 2>H owing to the lin e arized 

Einstein equation = kT^. Second, in contrast to the current paper. IHahn &: Paraniapd (120161 1 
did not single out the very important contribution to 6Tq, namely, the second term in the rhs 
of (j2.13p . which is directly proportional to <1> (see, however, t heir Appendix C, wher e the a uthors 
address this issue along with the connection to the approach of IChisari &: Zaldarriaeal (j201l|)). The 
mentioned term is absolute ly necessary for satisfying the perturbed energy conservation equation 
(jEingorn fc BrilenkovI 120151 ) and leads to the screening length A (j3.5p irrespective of the velocity- 
dependent contribution. 


3.3. Yukawa Interaction and Zero Average Values 


It is important to stress that, as a manifestation of the superposition principle, the second 
term in ()2.40p represents the sum of Yukawa potentials 




KC 


mr 


Svra r — r„ 


■exp(-g„) = - 


GNTUr. 


R — R„ 


■ exp 


|R — R,i 

A 


(3.10) 


coming from each single particle, with the same interaction radius A. Such a favourable situation 
is possible owing to the last term in the left-hand side (Ihs) of Eq. (I2.15p . whic h has been disre¬ 
garde d in (jEingorn &: Zhukll2012l l by mistake and erroneously compensated in (|Eingorn &: Zhuk 
2014 1 by inhomogeneous radiation of unknown nature. Actually, such radiation must not only 
possess negligible average energy density (requiring additional questionable reasoning), but also 
exchange the momentum with the nonrelativistic pressureless matter, despite the fact that no 
non-gravitational interaction between these two constituents has been assumed, and therefore the 
energy-momentum interchange is strictly forbidden. Here the mentioned unpardonable omission is 
rectified: the ill-starred term is reinstated, and there is no necessity in any additional interacting 
Universe components at all. 

The sum 4>n is certainly convergent at all points except at positions of the gravitating 

n 

masses, and computational obstacles do not come into existence. In particular, the order of adding 
terms corresponding to different particles is arbitrary and does not depend on their locations. On 
the contrary, there are certain obstacles when calculating the sum of Newtonian pot entials or thei r 
gradients. Let us address the well-known formulas (8.1) and (8.3) in the textbook bv IPeeblesI (|l98fll ) 
for the gravitational potential and the peculiar acceleration, respectively, derived in the Newtonian 
approximation (see above): 


4 > 


dr' 


/ P\r=r' P 


-V$ 


r — r 




r — r 


(3.11) 


up to space-independent factors being of no interest here. Substituting ()2.12p into the second 
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integral in (13.lip , one gets the formula (8.5) in (|Peebles 1980 1: 






V -Tr 


:(r-rri 


(3.12) 


According to [Peebles! (jl98d l. this sum is not well-defined and depends on the order of adding terms, 
and if one adds them in the order of increasing distances |r—r„| and assumes that the distribution of 
particles corresponds to a spatially homogeneous and isotropic random process with the correlation 
length being much less than the Hubbl e radius c /ff, t hen this sum converges. As regards the hrst 
integral in ()3.1ip . the argumentation by [Peebles] (|l980l l again relies on the random process assuring 
convergence; however, substitution of (I2.12p splits this integral into two divergent parts: the sum 
of an infinite number of the Newtonian potentials of the same sign (see also the paper by [Norton 
(j 19991 ) devoted to the related famous Neumann-Seeliger gravitational paradox) and the integral of 
the pure Newtonian kernel. 

Of course, the enumerated difficulties are absent when summing up the Yukawa-type potentials. 
In addition, it is interesting that in this case the particles’ distribution may be nonrandom and 
anisotropic. The lattice Universe model with the toroidal topo logy T xT xT represents a striking 
example. As explicitly demonstrated bv [Brilenkov et al.[ ([20151 1 . in the framework of this model the 
gravitational potential has no definite values on the straight lines joining identical point-like masses 
in neighbouring cells if the last term in the Ihs of Eq. (|2.15p is not taken into account. Evidently, 
the finite Yukawa interaction range A arising due to this term easily resolves this challenge as well 
as any similar ones related to the choice of periodic boundary conditions. Incidentally, if the space 
is supposed to have the usual, non-toroidal topology Rx Rx R, but the choice of periodic boundary 
conditions is made for N-body simulation purposes, then the dimensions of a cell should normally 
be greater than A, thereby weakening the undesirable impact of periodicity on simulation outputs. 


A noteworthy feature of the Yukawa potentials (I3.10p consists in assuring the zero average 
value of the scalar perturbation $ (|2.40p . Let us determine the average value of a single one of 
them: 


4 / 

V 


dripn = - 


KC rUn 
Sira V 


dr 


r - r. 


exp 


V 


a\r - Vr, 


KC rrin 


dvrA^ 


Sira V 


m„ 1 


where the comoving averaging volume V tends to infinity. Here the definition of A (13.5p has been 
used. Consequently, 

n ^ n 

since {1/V)^mn = p. Combining (|3.14p with the first term in (I2.40p . one immediately achieves 

n 

the desired result <I> = 0 (the third term in ()2.40p is apparently zero on average in view of the 
different directions of particle velocities, and the same applies to the vector perturbation B (j2.36l) : 
B = 0). This result means that the first-order backreaction effects are absent, as it certainly 
should be. Zero average values of the first-order cosmological perturbations are expected from 
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the very beginning, since these metric corrections are none other than linear de viations from the 
unper t urbed average v alues of the metric coefficients. Nevertheless, as shown bv lEingorn &: Zhuk 
(|2014l l: lEingorn et ahl (j2015l b there exists a concrete example of the mass distribution, which gives 
the nonzero average value of the gravitational potentia l determined by the standard prescription 
(I3.1ip . This problem is solved bv lEingorn et ahl (j2015l l through introducing manually the abrupt 
cutoff of the gravitational interaction range with the help of the Heaviside step function. One can 
see now that the same problem is strictly solved with the help of the hnite Yukawa range, and the 
potential remains smooth together with its gradient thanks to the smoothness of the exponential 
function. Obviously, the establi shed equality = 0 t akes place for an arbitrary mass distribution 
including that investigated by lEingorn et ahl ([201a). In addition, the well-grounded equalities 
6Tq = 0 and ST^ = 0 are valid as well, following from (I2.13P and (12.141) . respectively. 

Let us bring up and settle a related issue consisting in the following. One can easily prove 
that in the limiting case of the homogeneous mass distribution 4> = 0 at any point, as it certainly 
should be. For example, on the surface of a sphere of physical radius R the contributions from 
its inner and outer region s combined with 1/3 in (12.401) give zero (see, e.g., the expression (3.12) 
in ([Eingorn &: Zhukll201[ll l for the gravitational potential within a spherical shell of uniform mass 
density, the inner radius —)• 0 and the outer radius R 2 + 00 , which gives —1/3, with the 

exception of the zero mode, which should be dropped). This means that in the considered limiting 
case the equation of motion of a test cosmic body reads: 


R = -R, (3.15) 

a 

so the acceleration of the body is reasonably connected with the acceleration of the Universe 
expansion. At the same time, the described simple, but crucial test cannot be passed by Newtonian 
gravity. Indeed, in the framework of the Newtonian cosmological approximation the contribution 
from the outer region of the considered sphere is absent, while the contribution from its inner 
region generates an additional force in the rhs of Eq. (I3.15p . spoiling the established correspondence 
between the accelerations. This demonstrates once again the superiority of the formula (|2.40p for 
all scales. 


3.4. Transformation of spatial coordinates 


When writing down the perturbed metric (12.41) . the gauge choice is made in favour of the 
so-ca ll ed Poisson/longitud inal/conformal-Newtonian gauge, by analogy with lAdamek et ahl (|2013l . 


20141 1: iMilillo et ahl (|2015l l. However, it is common knowledge that there is no preferable coordinate 


system, so other gauges are admissible as well. The chosen gauge is characterized, in particular, by 
the coi ncidence of the found function <I> (|2.40p with the corresponding gauge-invariant Bardeen po¬ 
tential ( BardeenlllQSOl l. The introduced energy-momentum fluctuations 6T^ also coincide with the 
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corresponding gauge-independent quantities. For instance, let us verify that t he expressiori (12.131) 
for JTO remains unchanged for the analogue of the so-called A^-body gauge (jFidler et al.l l2015l i . 
This particular gauge features the unperturbed comoving volume giving a chance of eliminating 
the second term in the rhs of (j2.13p and, hence, of rehabilitating the Newtonian description. In 
this connection, it is necessary to show directly that this chance does not contradict the Yukawa 
screening of the gravitational interaction established in the Poisson gauge. For this purpose, let us 
rewrite the metric (j2.4p excepting the vector perturbation B: 


ds^ = 


(1 -|- 2^)dr]^ “ (1 “ 2^)5apdx°‘dx^ 


(3.16) 


where the scale factor a is a function of the conformal time r/ while the scalar perturbation (12.401) 
is a function of rj and comoving coordinates a = 1, 2, 3. The transformation of coordinates 


^ a dL 

V = T + A, x=e+ —, 


(3.17) 


where A and L are (first-order) functions of the new conformal time r and new comoving coordinates 
a = 1, 2,3, gives 


ds^ = 


(1 -F 24> -F 2A' -F 2'HA) dr^ -F 2 




- (1 - 24> + 2nA)6a^ + 2 


02 L 






drde 


(3.18) 


Here the prime denotes the derivative with respect to r; a and Ti depend on r while depends on 
(r, ^"). Fixing H = 0, one immediately comes to the opportune coincidence of the fluctuations of the 
mixed energy-momentum tensor components with the corre sponding gauge- inv ariant perturbations. 
Despite the fact that this choice differs from that made bv iFidler et al.l (120151 ) (where A^ 0), this 
does not affect the following main idea of the Y-body gauge. In accordance with the general 
definition (j2.1ip . in the new coordinates (t, ^") instead of (I2.13P one has 


2 — 2 

jpo = ^ (34> - AY) 


where = 5'^P 


92 


d^»d^P 


^ P^- P, Pi = '^ [c- fr 


(3.19) 


(3.20) 


Next, fixing A^L = 3<h, one may present the energy density fluctuation (I3.19|) conformably in the 
form 


ST^ = . 


(3.21) 


Thus, it may seem that proper use of gauge freedom ensures disappearance of the second term in 
the rhs of Eq. (I2.13|) . Nevertheless, the expressions (j2.13p and (I3.21|) for 5Tq are equal. In order 
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to prove this, let us use the fact that the perturbation 6p^ entering into (I3.2ip is not equal to the 
counterpart 6p entering into (12.131) . Indeed, the rest mass density p (I2.12p is connected with p^ 
(j3.20p by means of the relationship 

(3.22) 


P = 


-Pi 


1 + 

where the denominator represents the Jacobian det of the comoving coordinates trans¬ 

formation. Since p = -p + 6p and pj = p + Jpg, recalling that L is the first-order quantity, from 
(j3.22p one gets 

6p^ = 5p + -pA^L = 6p + 3p<I>. (3.23) 

Substitution of ()3.23p into (13.211) revives the gauge-independent perturbation (I2.13p . It is important 
to remember that positions of the gravitating masses are described by radius-vectors which certainly 
depend on the choice of comoving coordinates. For instance, apparently, ^ in the case of the 
nontrivial function L in ()3.17p . 


The initial displacement of particles proposed bv IChisari fc Zaldarriagal (120111 ) can be studied 
in the same vein. Restricting themselves to the linear relativistic perturbation theory for large 
enough scales where the failure of Newtonian dynamics is expected and striving for absorption 
of relativistic effects into the initial conditions for Newtonian simulation codes, the authors took 
advantage of the transformation of spatial coordinates 

d 


+ 5x 


a 

in’ 




(Jx“) = 3Ci, 


(3.24) 


where On stands for th e initial value of the so-called comoving curvature, or curvature perturbation 
variable (jPurreii 120081 ). 

+ (3,25) 


KpC^ 


Then substitution of (|3.19p . where now A^L is replaced by 3Cin, into (12.71) gives 


A$ - + n^) = ‘^ [5p^ + 3p ($ - Cin)]. 

la 


(3.26) 


Taking into account that the introduced comoving curvature does not evolve at large scales under 
consideration, one can replace Cin in (13.261) by C (|3.25p . and the subsequent cancellation of terms 
in the obtained equation reduces it to the following form: 


Ad> =- 5pc 

2a 


(3.27) 


Once again, as it follows from the first equality in ()3.23l) . 6p^ = 6p + 3pCm- Then Eq. (I3.27p is 
reduced to its original form before the transformation (j3.24p . in complete agreement with the gauge 
invariance of the Bardeen potential. 

Summarizing, there are two consistent options for cosmological simulations. O n the one hand, 
one can resort to the initial displacement of particles (jChisari fc Zaldarriagal 120111 ) or the A-body 
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gauge ( Fidler et alJboiS ) and reinterpret the large-scale Newtonian N'-body outputs as the relativis¬ 
tic ones. On the other hand, one can remain faithful to the Poisson gauge and c alculate the gravita¬ 
tional potential from the Hel mholtz equation, in harmony with the reasoning bv IHahn &: Paraniaoe 
(|201a i (see also the paper by iRampf &: Rigopoulosl (|2013l l where the Helmholtz equation links the 
potential to the density perturbation at scales comparable to the horizon). 


3.5. Nonzero Spatial Curvature and Screening of Gravity 


The promised generalization to both cases of curved spatial geometry can be made straightfor¬ 
wardly. For simplicity and illustration purposes, let us restrict ourselves to Eq. (I2.27P and rewrite it 
dropping the velocity contributions (i.e. the second term in the rhs) and taking into consideration 
the nonzero spatial curvature: 


A<^ + 



2>K,pc^ \ 
2a ) 




(3.28) 


where /C = -|-1 for the spherical (closed) space and /C = — 1 for the hyperb ol ic (open) space , and t he 


Laplace operator is redefined appropriately ('see Eingorn & Zhuk ( 

201 

2): I 

Bureazli et al. 

(2015)). 

This equation is equivalent to the equation (2.25) in ( 

Bureazli et al. 

201 

5) up to desig 

nations, 
ting the 

Hence, one can make use of its solutions derived bv Bureazli et al 

(2 

015) 

, simply adius 


notation. There seems no sense to reproduce these solutions here, but it should be emphasized that 
they are smooth at any point except at particle positions (where the Newtonian limits are reached) 
and characterized by zero average values, similarly to the flat space case /C = 0. 


One more important detail lies in the fact that the definition of A (|3.7p remains valid not only 
in the curved space case, but also in the presence of an arbitrary number of additional Universe 
components in the form of perfect fluids with constant or varying paramete rs in the equations of 
state like p = toe (e.g., radiation with the parameter 1/3), as one can prove (jEingorn fc Brilenkov 
20151 ). This hints at the universali t y of the presentation (13.7p . In particular, the gravitational 


potentials derived bv lBurgazli et al.l (1201.51 ) may be interpreted as valid for the Universe filled with 
quintessence with the parameter —1/3 in the presence of the cosmological constant as well as the 
nonrelativistic pressureless matter with negligible average energy density. 


Returning to the conventional cosmological model, from (|3.7I) one gets the dependence A ~ 
at the radiation-dominated stage of the Universe evolution. Since A may be associated with the 
homogeneity scale, as stated above, the asymptotic behavior A —>■ 0 when a —)• 0 supports the idea 
of the homogeneous Big Bang. 


Finally, it seems almost impossible to overcome the irresistible temptation of associating 
the Yukawa interaction range A with the graviton Compton wavelength h/{mgc), where h is 
the Planck constant and nig is the graviton mass, in the particle physics spirit. However, one 
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should act with proper circu mspection when discussing the massive graviton (se e reasoning by 


Faraoni &: CooperstockI (|l998l i as well as argumentation bv iGazeau &: Novellol (|201lh with respect to 


Minkowski and de Sitter spacetimes). It is remarkable that setting A equal to h/{mgc), h = hj{2iT), 
gives rUg = h/ (Ac) ss 1.7x 10“^^ eV today (when a = qq), and the ratio 1/A^ = m?gC^/f^ turns out to 
be numerically equal to 2A/3. And vice versa, if one does not initially resort to the known numerical 
valu e oi Qm and, hence, does not estimate A and rUg, the conjectural relationship rrigCp= 2A/3 
(see iHaranas fc Gkigkitzisl (120141 ) and Refs, therein) when a = uq may be rewritten with the help 
of (13.5p as 9 Gm = 4Ga, whence in the case of the negligible spatial curva ture {Q-m + Ga = 1) 
one gets VLm = 4/13 ~ 0.31, in solid agreement with lAde et al.l (j2014l . l2015l l. It is noteworthy as 
well that since A ~ (|3.5I) . one has rUg ~ so rUg ~ 1/t at the matter-dominated stage 

of the Universe evolu t ion (w hen a ~ This dependence on time agrees with that found by 


Haranas fc GkigkitzisI (120141 1. At the radiation-dominated stage A ~ a^, nig ~ a Thus, nig 
when a — )• -|-oo (A-domination prevents screening of gravity) and formally nig — -|-oo when a 


> 0 

0 . 


The established finite Yukawa range of the gravitational interaction may potentially pretend 
to play a key role in resolving the coincidence and cosmological constant problems as well as 
developing the holographic and inflationary scenarios. Glarihcation and rigorous substantiation of 
this role overstep the limits of the current paper. 


4. CONCLUSION 

The following main results have been obtained in the present paper in the framework of the 
concordance cosmological model: 

• the first-order scalar (I2.40p and vector ()2.36p cosmological perturbations, produced by inhomo¬ 
geneities in the discrete form of a system of separate point-like gravitating masses, are derived in the 
weak gravitational field limit without any supplementary approximations (no 1 /c series expansion, 
no “dictionaries”); 

• the obtained explicit analytical expressions (12.361) and (12.401) for the metric corrections are valid 
at all (sub-horizon and super-horizon) scales, converge at all points except at locations of the 
sources (where the appropriate Newtonian limits are reached), and average to zero (no hrst-order 
backreaction effects); 

• both the Minkowski background limit (see (|3.2p and (13.31) 1 and the Newtonian cosmological 
approximation (see (I3.4p ). which is widely used in modern N-body simulations, represent particular 
limiting cases of the constructed scheme and serve as its corroboration; 

• the velocity-independent part of the scalar perturbation (I2.40p contains a sum of Yukawa poten- 
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tials produced by inhomogeneities with the same finite time-dependent Yukawa interaction range 
(I3.5p . which may be connected with the scale of homogeneity, thereby explaining the existence of 
the largest cosmic structures; 

• the general Yukawa range definition (13.71) is given for various extensions of the ACDM model 
(nonzero spatial curvature, additional perfect fluids), and advantages of the established gravity 
screening are briefly discussed. 


Based on the obtained results, it should be not too difficult to construct similarly an appropriate 
scheme for the second-order cosmological perturbations including the tensor ones. Accomplishment 
of this quite possible technical mission would predict, in particular, the backreaction effects. It is 
expected that the second-order metric corrections will be much smaller than the first-order ones 
at arbitrary scales. Besides, the direct generalization of the elaborated approach to the case of al¬ 
ternative (nonconventional) cosmological models, for example, those replacing the A-term by some 
other dynamical physical substance, serving as dark energy and also fitting all data, is straightfor¬ 
ward and can be made with hardly any trouble. Then, simulating nonlinear dynamics at arbitrary 
scales, predicting formation and evolution of large cosmic structures, determining the influence of 
metric corrections on propagation of photons through the simulation volume, etc, one can actually 
probe cosmology and potentially distinguish among different competing representations of the dark 
sector. Of course, extra effort and care are required for constituting a link between physical quan¬ 
tities extracted from relativistic simulation s and observables measu r ed in galaxy surveys such as 
redshifts and positions in the sky (see, e.g.. iBonvin &: Durreii (j201lh : lYoo &: Zaldarriaeal (|2014l B. 


Thus, the developed cosmological perturbations the ory covering the whole space, in combi¬ 
nation with such future high-precision surveys as Euclid (jScaramella et al.l 120151 ). approaching the 
Hubble horizon scale, may essentially deepen our knowledge about the amazing world we live in. 
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